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Abstract. Convergence problems in coupled-cluster iterations are discussed, and 
a new iteration scheme is proposed. Whereas the Jacobi method inverts only the 
diagonal part of the large matrix of equation coefficients, we invert a matrix which also 
_J | includes a relatively small number of off-diagonal coefficients, selected according to the 

O ' excitation amplitudes undergoing the largest change in the coupled cluster iteration. 

c/3 , A test case shows that the new IPM (inversion of partial matrix) method gives much 

£^~* better convergence than the straightforward Jacobi-type scheme or such well-known 

convergence aids as the reduced linear equations or direct inversion in iterative subspace 
methods. 






> 

o 
\o 
o 

On 
O 
O 
O 

O 



X 



PACS numbers: 31.15.Dv, 31.15.-p, 31.25.-v, 31.15.Ar 



Submitted to: J. Phys. B: At. Mol. Opt. Phys. 



E-mail for correspondence: Mosyagin@lnpi.spb.su; http://www.qchem.pnpi.spb.rvj 



Letter to the Editor 2 

The coupled cluster (CC) method is widely used in electronic structure calculations. 
The CC theory has been described in many reviews (see, e.g., |], 0, [3|, |J |5]), and will 
not be presented here. The basic equation for the CC method is the Bloch equation 

nun = un, (i) 

where H is the Hamiltonian and f2 is the wave operator. The resulting equations have 
the general algebraic form 

TV 

A i + J2B(t) i3 t,=0 1 2 = 1,2,. ..,N, (2) 

i=i 

where tj are the cluster or excitation amplitudes to be determined, iV is the number of 

the unknown amplitudes, A is a vector and B(t) is a square matrix which in general 

depends upon t. For simplicity, we consider the case when B does not depend upon t, 

N 

A i + Y / B lJ t J = 0, i = l,2,. ..,N. (3) 

3=1 

The generalization for the case of B(t) is straightforward (and implemented in the 
relativistic CC code employed for the test examples below). The direct solution of 
equations (|3|) using the Gauss elimination method is feasible only for systems with a few 
thousand cluster amplitudes at most, whereas problems encountered in our relativistic 
CC may involve millions of such amplitudes. A Jacobi-type iterative method is usually 
applied to solve these equations. Using the fact that B is normally a diagonally dominant 
matrix, the method involves direct inversion of the diagonal part D of B. The system (|3|) 
is rewritten in the form 

N 

U = -(D- 1 )^ + J2(B - D^tj], i = l,2,...,N (4) 

3=1 

and is solved iteratively. 

The coupled cluster calculations are often beset by convergence difficulties. This is 
particularly true for multireference CC methods, such as the Fock-space approach [|2|. ||]. 
Several methods for improving convergence have been proposed; the most commonly 
used are the reduced linear equations (RLE) || and direct inversion in the iterative 
subspace (DIIS) [0, |J approaches. These help in some, but not all, cases. Most severe 
convergence problems may be traced to the existence of intruder states. While increasing 
the model (or P) space improves the quality of the calculation by including a larger part 
of the correlation, it also increases the probability of encountering intruder states and 
getting no valid results at all. New methods for improving convergence are therefore 
highly desirable. One such method is presented in this Letter. 

The problem may be illustrated by an example taken from recent work ||, where 
ground and excited state energies of Hg and its ions were calculated by the relativistic 
coupled cluster method. The 5<i 10 ground state of the Hg 2+ ion served as the reference 
state, and the Fock-space CC scheme was 

Hg 2+ [(0)sector] -► Hg+[(l)sector] -► Hg[(2)sector], (5) 
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with electrons added in the 6s and 6p orbitals, designated as valence particles. While 
the calculations in |J were relativistic, the nonrelativistic notation will be employed 
for brevity. The model space in the (1) sector, with one valence particle, consisted 
of determinants with 5d 10 6s x and 5<i 10 6p 1 configurations. Adding the 7s, 7p, and 6d 
orbitals to the list of valence particles would yield more state energies in the (1) sector, 
as well as better description of the 6s 1 6p 1 states |J in the (2) sector, corresponding to 
neutral Hg. Unfortunately, adding the 5d 10 7s 1 , 5d 10 7p 1 , and 5d 10 6d 1 configurations 
to the model space leads to divergence of the CC iterations (U) in the (1) sector. 
Analysis shows that the divergence is caused by the 5d 9 6s 1 6p 1 , 5d 9 6s 2 and other intruder 
states from the complementary Q space, which are close in energy to certain P states 
(5d 10 7p 1 , 5d w 6d 1 , and others). The diagonal elements B it of the matrix B correspond 
to differences between the total energies of the P and Q determinants connected by the 
ti excitations. Some of these elements will be very small in this case, leading to large 
elements in D^ 1 . Small changes in t amplitudes on the right hand side of equations (^) 
will therefore cause large changes in the amplitudes on the left hand side, leading to 
divergence. 

We propose to overcome this problem by replacing the D matrix by D' which 
includes, in addition to the diagonal elements of B, those nondiagonal B elements which 
are large in comparison with corresponding diagonal elements. The calculation of all B 
matrix elements is impractical, and a selection procedure for nondiagonal elements to 
be included in D' is described below. This new matrix is constructed so that its matrix 
elements, D\ , are equal to or approximate the Bij matrix elements both for i = j and 
for i,j G /, where I is some small subset of the amplitudes. The other nondiagonal D\, 
matrix elements (i (jL I or j : G" I) are set to zero. The method involves the inversion of 
the partial matrix (IPM) D'. A modified form of the system of equations (HI), 

JV N 

U = - Y,{D'- x ) ik {A k + £(£ - D') kj t 3 }, i = 1, 2, . . . , N (6) 

fc=l j=l 

is obtained and solved iteratively. Equations (H) can be divided into two sets, 

N 

U = - ^(D'^hklAk + E( S " D )kjtj ~ £(£' - D) kj t 3 ], i G /, (7) 
fee/ j=i j<=i 

N 

U = -(D'- 1 )^ + J2(B - D^jtj - (D' - DhU], % £ /, (8) 

3=1 

where the second part is similar to equations (f|). 

The size of the subset I must be kept small, so that the calculation (the most time- 
consuming step), storage and manipulation of the non-zero off-diagonal D' elements 
remains feasible. Careful selection of the amplitudes to be included in I is therefore 
of paramount importance. The algorithm followed here starts with calculating the ti 
amplitudes by the standard iteration scheme (01). The amplitudes which have undergone 
the largest changes are included in /, the corresponding D'^ matrix elements are 
evaluated, and the ti amplitudes in I are recalculated by equations (0). The dimension 
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M of the I subset was kept at 1000, which makes the calculation and manipulation of 
D' feasible. Optimal algorithms for determining M, selecting excitations to be included 
in /, and calculating the D' matrix will be studied in the future. 

It should be noted that the system (^) is equivalent to the standard equations (f|) 
in the limit M = 0; in the limit M = N, scheme ([]) converges in one iteration, if one 
takes D[j = Bij, 

N 

t i = -Y,(B- 1 ) ik A k , i = l,2,...,N. (9) 

fc=i 
Formally, one can always achieve convergence of the iterations (|]) by increasing M. 
The IPM method proposed here may be combined with other procedures for accelerating 
convergence, such as the reduced linear equations || and direct inversion in the iterative 
subspace [0, §] methods. This has not been done in the present application, and will 
be tried in the future. It should be mentioned that the identification of the resulting 
high-lying levels may require careful analysis of the t amplitudes, particularly if some 
of the latter are large, indicating large contributions of Q configurations. Finally, it 
should be noted that the IPM scheme described above may be regarded as adopting the 
Gershgorn-Shavitt A\. perturbation theory approach rather than that of A |T0] . 



The different iteration schemes were tested for the 33-electron relativistic Fock- 
space CC calculation with single and double cluster amplitudes of Hg + levels in the 
{spdfg) basis from || in the framework of the Dirac-Coulomb Hamiltonian. Two 
model spaces were used, one consisting of determinants with 5rf 10 6s 1 and 5c? 10 6p 1 
configurations, the other including in addition the 5<i 10 7s 1 , bd 10 7p l , and bd w Qd l 
configurations. All iterations involved 1:1 damping (the input amplitudes for iteration 
n + 1 were taken as the average of input and output amplitudes of iteration n). The 
IPM scheme is compared with the standard scheme (§) and with the RLE || and DIIS 
0, H methods in tables [I] and |[ The RLE and DIIS methods used the output of 
the last five iterations to form the new input vector. All methods led to convergence 
for the small model space (table [I]). The RLE, DIIS, and IPM schemes were about 
equally effective in reducing the number of iterations required. The large model space 
(table ^) shows markedly different behavior for the different methods. Straightforward 
iteration by the Jacobi-type method blows up almost immediately; the large excitation 
amplitudes may be traced to the intruder states mentioned above. The RLE and DIIS 
schemes exhibit better behavior, but could not achieve convergence even after several 
hundred iterations. Only the IPM approach proposed in this Letter led to convergence 
(in the 29th iteration), showing the potential of the method. 
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Tables and table captions 



Table 1. 

N 



The largest change in the single and double cluster amplitudes 



(max|td™ — 4 I) a * iteration n. The changes are obtained by equations (0) in 

i— 1 LJ 

the RCC calculations with the Jacobi-type, RLE, DIIS and IPM iteration schemes. 
The model space consists of determinants with 5d 10 6s 1 and 5d 10 6p 1 configurations. 
The convergence threshold is 10~ 6 . 



Iteration 


Jacobi 


RLE 


DIIS 


IPM 





1.62- lO" 1 


1.62- IO" 1 


1.62- IO" 1 


1.62- IO" 1 


3 


2.99 


io- 2 


2.99- IO" 2 


2.59 


io- 2 


1.30- IO" 2 


6 


1.39 


io- 2 


1.41 • IO" 3 


1.10 


io- 3 


1.38 -IO- 3 


9 


6.65 


10~ 3 


6.81 • IO" 4 


1.21 


io- 4 


2.01 • IO" 4 


12 


3.19 


10~ 3 


1.69- IO" 5 


1.12 


io- 5 


3.60- IO" 5 


15 


1.54 


10~ 3 


8.42 ■ IO" 6 


5.15 


io- 6 


7.01 • IO" 6 


18 


7.53 


io- 4 


9.47- IO" 7 


2.43 


io- 6 


1.43- IO" 6 


21 


3.72 


io- 4 


convergence 


1.15 


1Q- 6 


convergence 


24 


1.86 


io- 4 




convergence 




27 


9.46 


io- 5 








30 


4.87 


io- 5 








33 


2.53 


io- 5 








36 


1.34 


io- 5 








39 


7.11 


io- 6 








42 


3.82 


IO- 6 








45 


2.07 


io- 6 








48 


1.13 


io- 6 










convergence 
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Table 2. Same as Table 1, except that the model space is larger, consisting of 
determinants with 5d 10 6s 1 , 5d w 6p 1 , 5d 10 7s 1 , 5d 10 7p 1 , and 5d w 6d 1 configurations. 



Iteration 


Jacobi 


RLE 


DIIS 


IPM 





2.14- lO" 1 


2.14- IO" 1 


2.14 


io- 1 


2.14- IO" 1 


3 


7.54 -10- 1 


7.54 • 10- 1 


6.93 


io- 1 


9.38 


io- 2 


6 


2.61 


5.18 - 10" 1 


1.22 


io- 1 


1.21 


io- 2 


9 


9.31 


1.82 


5.75 


io- 2 


1.25 


io- 3 


12 


3.42 • 10 1 


2.73 • lO" 1 


3.20 


io- 2 


3.19 


io- 4 


15 


1.24 


10 2 


5.65 • 10- 1 


2.90 


io- 2 


7.39 


io- 5 


18 


4.45 


10 2 


1.80 -10 1 


2.87 


io- 2 


1.61 


IO- 5 


21 


1.54 


10 3 


2.77 -HT 1 


3.06 


io- 2 


1.04 


io- 5 


24 


5.00 


10 3 


1.01 


2.55 


io- 2 


9.55 


io- 6 


27 


1.49 


10 4 


1.46 • IO" 1 


2.19 


io- 2 


3.64 


1Q- 6 


30 


3.93 


10 4 


5.51 • 10- 1 


2.23 


io- 2 


convergence 


33 


9.08 


10 4 


2.40 • IO" 1 


1.93 


io- 2 




36 


1.88 


10 5 


1.73 • 10- 1 


1.99 


io- 2 




39 


3.61 


10 5 


4.91 • lO" 1 


1.36 


io- 2 




42 


6.70 


10 5 


2.95 • 10- 1 


1.62 


io- 2 




45 


1.22 


10 6 


6.04- IO" 1 


1.44 


io- 2 




48 


2.22 


10 6 


5.00 • 10- 1 


1.45 


io- 2 






divergence 


no convergence 


no com 


mergence 





